Classical simulation of quantum algorithms using the tensor product representation 



IT) . 

o ■ 
o ' 

<n : 

>v 
<3 



m 

(N 
> 

in 
o 

(N 



■ 

o ■ 

Ok 
^— > ■ 

G ■ 

3 : 



A. Kawaguchi 1 ' 3 , K. Shimizu 1 , Y. Tokura 1 ' 3 and N. Imoto 1 ' 2,3 ' 4 
1 NTT Basic Research Laboratories, NTT Corporation, 
3-1 Morinosato-Wakamiya, Atsugi, Kanagawa 243-0198 
2 Department of Materials Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, 
3 SORST Research Team for Career Electronics, JST 
4 CREST Research Team for Photonic Information Processing, JST 
(Dated: February 1, 2008) 



Japan 




Using the tensor product representation in the density matrix renormalization group, we show 
that a quantum circuit of Grover's algorithm, which has one-qubit unitary gates, generalized Toffoli 
gates, and projective measurements, can be efficiently simulated by a classical computer. It is 
possible to simulate quantum circuits with several ten qubits. 



Quantum computation is currently attracting much in- 
terest Shor's factorization algorithm [3j, which 
finds the prime factorization of an ro-bit integer in the 
polynomial time of n, proved the power of quantum com- 
putation . Grover's search algorithm which has com- 
plexity 0{\/N) instead of 0(N) classical steps, where 
N is the database size, provides proof that quantum 
computers are faster than classical ones for unstructured 
database searches. Moreover, it is expected to have ap- 
plications to other important NP problems. 

Actual physical systems are affected by various deco- 
herence and operational errors, which make it difficult to 
analyze the properties of quantum circuits. Therefore, 
it is very important to simulate a quantum computer 
realized with physical systems by a classical computer. 
However, simulating large quantum circuits with a classi- 
cal computer requires enormous computational power be- 
cause the number of quantum states increases exponen- 
tially as a function of the number of qubits. Recently, Vi- 
dal presented an efficient method for simulating quantum 
computations with one-qubit unitary gates and neighbor- 
ing two-qubit unitary gates Q. Before that, White had 
introduced a numerical technique to study the proper- 
ties of relatively large-scale one-dimensional (ID) quan- 
tum systems, which is called the density matrix renor- 
malization group (DMRG) 0. Vidal's idea and DMRG 
strongly overlap. As an interesting combination of the 
above two methods, a time-dependent DMRG method 
has been presented H,@. 

In this paper, we present a new simulation method 
that uses the graphical representation of the tensor prod- 
uct state 0. With this method, we can take advantage 
of arbitrary one- and two-qubit operators, which include 
projection operators, and can execute the projective mea- 
surement in an efficient way. In addition, generalized 
Toffoli gates (C fc -NOT) can be used without decomposi- 
tion in terms of elementary gates. As an example of the 
simulation, we will show that Grover's algorithm can be 
efficiently executed by a classical computer. Here, to im- 
plement Grover's algorithm concretely, we assume that 
the Oracle consists of generalized Toffoli gates. Using 
this simulation method, it becomes possible to simulate 



FIG. 1: Graphical representation of the tensor product state 
for L = 6 in Eq. JTJ. 



various quantum circuits with several ten qubits. 

According to the tensor product decomposition, any 
wave function of pure state \I/o with L qubits can be 
decomposed to L tensors. 

m 2 mi 

*o(*i,*2, ■■■ ,ih) = ^2 ■■■ ^2 

xi = \ x L = l 

Ui(ii,xi,X2)U 2 (i2,X2,x 3 ) ■ ■ ■ U L {i L ,x L ,x L+ i), (1) 

where i a represents the a-th qubit. The right side con- 
sists of the product of L tensors, and x a {= 1, • • • m a ) are 
new variables related to the entanglement between neigh- 
boring qubits 0- Here, X\ = x^+i = 1 (mi = m^+i = 1) 
for the pure state. This representation contains 0(m 2 L) 
components, where m is the maximum value of m a . 

In Fig. ^ we show the graphic representation of Eq. 
with L = 6 6J. In this graph, we represent tensors 
U a as triangles, and each vertex corresponds to i a (cir- 
cles) and x a (squares). We use black marks in the same 
manner as Einstein's summation convention. The repre- 
sentation implies that a summation is taken over the cor- 
responding index when an index is repeated once in the 
same term. Here, all m a become m a = 1 for the product 
states. Therefore, tensors U a (i a ,x a ,x a +i) are reduced 
to vectors U a (i a ). For example, in the case of a superpo- 
sition state, \l1 Ll2 Y?i=i l*)> then U » = (1/V2, 1/V2) T . 
Also, when a state is 1 00 • • • 0), then U a — (1, 0) T . In this 
paper, we use the 1 00 • • - 0) state as the initial state to 
simulate quantum circuits. 

As a simple example, we explain a calculation process 
for the inner product. Fig. 121 is the graphical represen- 
tation of the inner product (^ol^o) f° r L — 6, where 
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FIG. 2: Graphical representation of the inner product 
(*o|*o). 




FIG. 3: Graphical representation of pi(ii,Xi\ii,x'i) and 
P2(k 2 ,y 2 \k'2,y' 2 ). 

dotted lines imply that a summation is taken over the 
connected variable. We calculate (^ol^o) as follows: 

X G {x G ,x' & ) = ^2 U e (i 6 ,xe,xr)Ug(i 6 ,x' 6 ,xr) (2) 

ie,x 7 

Y 5 (i 5 ,x 5 ,x' 6 ) = ^2 X 6 (x 6 , x' 6 )U 5 (i 5 , x 5 , x 6 ) (3) 

xe 

X 5 (x 5 ,x' 5 ) = ^2 Y 5{k,X5,x' 6 )Ug(i 5 ,x' 5 ,x' 6 ) (4) 
X 1 (x 1 ,x' 1 ) = ^2 Yiih^t^^Uiih^'^x^). (5) 

21,^2 

Here, Xi(xi,x'i) has the scalar quantity because x\ = 
x{ = 1; therefore, (^ol^o) = Xi(l,l). The calculation 
of (^ol^o) f° r L qubits requires 0(Lm 3 ) multiplications. 

To simulate a quantum circuit one step further, we 
have to obtain the tensor product decomposition of the 
next state I*!) = O p \^> ), 

rn' 2 m' L 

*i(«x,i2, •••,*!,) = 2J *" zJ 

V2 = l VL = l 

Vi(ii,yi,y2)V 2 (i2,y2,y3) ■ ■ ■ V L (i L , y L , y L +i), (6) 



where y± — yL+i = 1- Here, we assume L = 6 and the 
operator is O p = 2 ^{k 2 , ki\i 2 , ii)- First, we consider the 
reduced density matrix p\ (%\ , x\ \i\ , x[ ) in Fig. IHJa) to get 
Vi and m! 2 . The order of the calculation becomes slightly 
complicated, but the calculation can be performed simi- 
larly to that for the inner product. Then, we diagonalizc 
the 2x2 density matrix pi(ix,xi\i' l ,x' x ), 

^2 v *( i uXi,y2)Pi(ii,x 1 \i' v x' 1 )V 1 (i' 1 ,x[,y 2 ) 

i\ ,x\,i'i ,x' r 

= ^y2,y' 2 ^V2 ■ CO 

Here, all eigenvalues \ V2 have positive values in the de- 
creasing order ( Ai > A 2 ). Then, we can rewrite 
Vi(ii,Xi,y 2 ) = Vi{ii,yi 1 y 2 ) since x\ = yi = 1 (mi = 
m'i = 1). Furthermore, we choose m' 2 as m' 2 — 2 when 
A 2 ^ 0, and m' 2 = 1 when A 2 = 0. When DMRG is used 
in the analysis of ID quantum lattice systems, such as 
the Heisenberg model, it is preferable to cut m' a that the 
(m a + l)-th eigenvalue of a density matrix becomes small 
enough. This corresponds to an irreversible compression 
that disposes of the unimportant information, and such 
compression is very effective. However, because quan- 
tum calculations, such as Grover's and Shor's algorithms, 
tend to treat discrete wave functions, important informa- 
tion is left out if we make m a too small. Therefore, we 
determine m' a such that the eigenvalue of a density ma- 
trix becomes in this paper. Namely, m' a corresponds 
to the Schmidt rank of the density matrix p a -i and the 
compression corresponds to the reversible compression 
that the original information is fully maintained. Next, 
we use V\ obtained previously and calculate the density 
matrix P2{k2,y2\k 2l y' 2 ) as Fig. CJb). Here, the qubit 
variable changes from i 2 to £; 2 because of 2 ,4- Then, 
V 2 and m' 3 can be obtained by diagonalizing p 2 in the 
same manner as obtaining V\ and m' 2 . By repeating the 
above-mentioned operations, we can get all V a and m' a . 

Next, we evaluate the computational complexity for 
one-qubit and two-qubit operations and show that the 
complexity is 0(Lm 3 ) for a projection operator. The 
calculation for obtaining p\ requires 0(Lm 3 ) multipli- 
cations, and 0(m 3 ) multiplications are necessary in or- 
der to diagonalize it. In the calculation for preparing 
p a (a > 1), some parts (ex. dashed boxes in Fig. 
|3J) are already obtained through the previous calcula- 
tion. This reduces the calculation for preparing p a to 
0(m 3 ) multiplications. Therefore, the total number of 
multiplications is 0(Lm 3 ) for getting all V a and m' a . 
Note that, for a one-qubit unitary Oj and a neighbor- 
ing two-qubit unitary operator Ojj+i, the calculation 
complexities can be reduced to 0(m 2 ) and 0(m 3 ), re- 
spectively 0] ■ In the case of the one-qubit unitary oper- 
ator Oj, V 3 {k 3l x 3l x J+1 ) = J2i 3 O j {k J1 i j )U l {i j ,x 3l x J+1 ) 
and V a = U a for a ^ j. Also, for the unitary oper- 
ator Ojj+i, the computational complexity for obtain- 
ing the next state becomes 0(m 3 ) by using V a = U a 
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(b) 



FIG. 4: Full Grover's algorithm (a) and Grover operator (b) 
for n — 4 with the number of solutions t — 2. In this example, 
correct solutions are 1 1011) and 1 1 100) . 



( a 7^ j,j + 1 ). However, for the projection operator, 
all U a may change due to the collapse of the wave func- 
tion. Thus, the computational complexity is 0(Lm 3 ) for 
a projection operator. We will explain the method of 
projective measurement later. 

Furthermore, we can use the C fc -NOT gate with k 
controlled qubits and one target qubit without decom- 
posing in terms of elementary gates. Though the G k - 
NOT gate can be decomposed into single qubit unitary 
gates and C-NOT gates 0, the C fc -NOT gate is very 
useful for constructing a quantum circuit. As a sim- 
ple example, we show the case of a C 2 -NOT operation 
for L = 4, which has two controlled qubits [i% and 1%) 
and one target qubit (ii). In this case, all elements 
of \&o(ii, *2, «3, 1) are replaced with $o(ii, *2, is, 0) when 
both i\ and %2 are |1). Therefore, we can represent the 
next state 1%, is, i^) as follows. 



-* (l,l,i3,i4) + *o(l,l,i3,i4). 



(8) 



Here, 14 means the reverse of 14. The i\ and ii in the 
second and third term are fixed to 1. We can make the 
reduced density matrixes using the tensor product rep- 
resentation and obtain V a and m' a one after another. In 
this way, we can use C fe -NOT with 0(Lm 3 ) multiplica- 
tions. Namely, we can use almost all operations appear- 
ing in quantum circuits with less than 0(Lm 3 ) multipli- 
cations. 

Grover's search algorithm involves repeated applica- 
tions of a quantum subroutine, which is called the Grover 
operator G. In Fig. 21 we illustrate a configura- 



§ 600 




FIG. 5: Execution time of a classical computer for Grover's 
algorithm up to n = 35. The inset shows the amplitude of 
a correct solution for n = 8. Filled circles and open squares 
indicate t — 1 and t = 2, respectively. The execution time 
indicates the computational time from k = to the first max- 
imum, k ~ 7r v/AT/i/4, in the inset. 



tion example of n — 4 with the number of solutions 
t = 2. Here, H and X represent the Hadamard gate, 

H=-^=(\ V and the NOT gate, X = (® ^ 

In this study, we assume that the Oracle [shaded part in 
Fig. E[b)] consists of t C"-NOT gates. 

When we apply k times of G to a superposition state, 
the state is represented as 



G k \^) = cos 



where 9 — 2 arcsin 



2k + 1 



2k + 1 



l/3>(9) 

(10) 
(11) 



Here, indicates a sum over all x that are solutions, 
and J2 X indicates a sum over all x that are not solutions. 
In the inset of Fig. we show the amplitudes of a correct 
state for t = 1,2 and n = 8. For convenience, we have 
fixed one of the solutions in |11 • • • 1) and have plotted 
the amplitude of 1 11 ■ ■ ■ 1) . Solid lines indicate analytical 



- ji. These results confirm that tensor 



lines, 

product representation can simulate Grover's algorithm 
precisely. 

A correct solution is provided by measuring qubits 
when k ~ tt^J N/t/4. In our numerical simulation, a 
projective measurement is performed as follows. First 
we calculate an expectation value (\E , |ni|^/) of i-th qubit 
and generate a random number r (0 < r < 1) to consider 
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the probabilistic nature of quantum measurement, where 
^' = f 1 ) ' ^ e a PP^ ^° ^ e * 1 u bit 

/P = |0)(0| if^lftil*) < r 

\P 1 = |1)(1| if(*|ni|*) > r. 1 j 

In Fig. we show the execution time of the quan- 
tum circuit illustrated in Fig 01 up to n = 35 using 
an Intel Xeon 2.4GHz processor. What has to be no- 
ticed is that the computation time increases when the 
number of answers increases. This behavior is com- 
pletely opposite to that for an ordinary search problem. 
Here, we discuss the computational complexity of a clas- 
sical computer. The circuit possesses 0(\/ N/t log 2 N) 
one-qubit unitary gates, 0(ty/ N/t) C fc -NOT gates, and 
log 2 N projective measurements. The calculation com- 
plexities per one gate are 0(m 2 ), 0(m 3 \og 2 N), and 
0(m 3 log 2 Y), respectively. Because the total computa- 
tional complexity is mainly dependent on C fc -NOT gates, 
it becomes 0(VtNm 3 log 2 N). We performed the fitting 
of Ay/N (log N) with a fitting parameter A (solid lines 
in Fig. EJ, and these fitting lines are consistent with the 
simulation results. Therefore, we see that m is not de- 
pendent on the qubit number n(= log 2 N). Actually, we 
have found that m has the relationship m < t + 1. The 
total computational complexity for finding one correct 
state can be expressed by 0(\/Nt 3 - 5 log 2 N). Inciden- 
tally, the required memory capacity is 0(t 2 log 2 TV). If t 
is much smaller than N, this computational complexity 
is approximately equal to that for a quantum computer. 
In this way, we have found that the tensor product repre- 
sentation is very useful in simulating Grover's algorithm. 

In summary, we have demonstrated the classical sim- 
ulation of a quantum circuit by using the graphical rep- 
resentation of the tensor product state. This method 
can take advantage of almost all gate operations, which 
include projection operators and the generalized Toffoli 
gates, with less than 0(Lm 3 ) multiplications. As a result 
of the simulation of Grover's algorithm, we have found 
that m, which determines the computational time and 
the memory size in a classical computer, is related to the 



number of solutions. This simulation method enables 
us to simulate various quantum circuits with several ten 
qubits. Therefore, we expect that this method will be 
helpful in analyzing quantum circuits with decoherence 
and operational errors. In particular, since this simu- 
lation method is quite effective for a one-dimensional 
quantum system, it will be useful for analyzing several 
important experimental approaches, which include a one- 
dimensional chain of quantum dots [lfj and the all-silicon 
quantum computer [111 [l2j • 
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computation. We also thank K. Okunishi and T. Nishino 
for discussions about the graphical representation of ten- 
sor product state. 
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